/******************************************************
Compile : gcc -O3 -o inferred1 inferred1.c sim.c solver.c utils.c -lm
Usage : ./inferred1 input-1
******************************************************/

#include "solver.h"
#include <stdio.h>
#include <string.h>

void print_pathway (int N, int L, int pathway[L][N])
{
	int i, j;

	printf ("c Candidate Pathway of Cell Cycle\n");
	printf ("%d %d\n", N, L);
	for (i = 0; i < L; i++)
	{
		for (j = 0; j < N; j++)
		{
			if (pathway[i][j] == -1)
				printf ("* ");
			else
				printf ("%d ", pathway[i][j]);
		}
		printf ("\n");
	}
}

/* restriciton of the turnover */
int network_inferred2 (int N, int L, int R[N + 1], int C[N][L - 1][R[N] + 1], int S[L][N], solver *s, int base)
{
	int i, j;

	veci lits;
	veci_new (&lits);
	lit *begin;
	int var;

	int T[L - 2]; /* turnover variable */

	/* restriction of the maximum R[i] turnovers for one nodees */
	for (i = 0; i < N; i++)
	{
		for (j = 0; j < L - 2; j++)
		{
			T[j] = base + j;
		}

		for (j = 0; j < L - 2; j++)
		{
			setup_turnover (S[j][i], S[j + 1][i], T[j], s);
		}

		base = counter_positive1 (L - 2, T, R[i], R[N], C[i], s, base + L - 2);
	}

	veci_delete (&lits);

	return base;
}

/* restriction for order */
int network_inferred3 (int N, int L, int MAX_R, int C[N][L - 1][MAX_R + 1], int M, int R2[M][4], solver *s, int base)
{
	int i, j;

	int X[L - 2], Y[L - 2];

	for (i = 0; i < M; i++)
	{
		for (j = 1; j < L - 1; j++)
		{
			X[j - 1] = C[R2[i][0]][j][R2[i][1]];
			Y[j - 1] = C[R2[i][2]][j][R2[i][3]];
		}
		base = check_order (L - 2, X, Y, s, base);
	}

	return base;
}

int main (int argc, char *argv[])
{
	FILE *fp;
	char buf[256];
	int i, j, m;
	int N, L, M;

	fp = fopen (argv[1], "r");
	fscanf (fp, "%d %d", &N, &L);

	int pathway[L][N], R1[N + 1];
	int S[L][N]; /* state variable */
	int g[N][N], r[N][N]; /* network variable ; g[i][j]/r[i][j] edge from j to i */

	for (i = 0; i < L; i++)
	{
		for (j = 0; j < N; j++)
		{
			S[i][j] = i * N + j;
			pathway[i][j] = -1;
		}
	}

	for (i = 0; i < N; i++)
	{
		for (j = 0; j < N; j++)
		{
			g[i][j] = L * N + i * N + j;
			r[i][j] = L * N + N * N + i * N + j;
		}
	}

	for (i = 0; i < N; i++)
	{
		fscanf (fp, "%d", &(pathway[0][i]));
	}

	R1[N] = 0;
	for (i = 0; i < N; i++)
	{
		fscanf (fp, "%d", &(R1[i]));
		if (R1[N] < R1[i])
			R1[N] = R1[i];
	}

	for (i = 0; i < N; i++)
	{
		if (R1[i] % 2 == 0)
			pathway[L - 1][i] = pathway[0][i];
		else
			pathway[L - 1][i] = 1 - pathway[0][i];

		pathway[L - 2][i] = pathway[L - 1][i];
	}

	fscanf (fp, "%d", &M);
	
	int R2[M][4];

	for (i = 0; i < M; i++)
	{
		fscanf (fp, "%d %d %d %d", &(R2[i][0]), &(R2[i][1]), &(R2[i][2]), &(R2[i][3]));
	}

	fclose (fp);

	solver *s;
	bool st;

	s = solver_new ();

	m = network_inferred1 (N, L, pathway, S, g, r, s);

	int C[N][L - 1][R1[N] + 1];

	m = network_inferred2 (N, L, R1, C, S, s, m);
	m = network_inferred3 (N, L, R1[N], C, M, R2, s, m);

	st = solver_simplify (s);
	if (st == true)
		st = solver_solve (s,0,0);

	if (st == true)
	{
		for (i = 1; i < L - 2; i++)
		{
			for (j = 0; j < N; j++)
			{
				if (s->model.ptr[S[i][j]] == l_True)
					pathway[i][j] = 1;
				else
					pathway[i][j] = 0;
			}
		}
		print_pathway (N, L, pathway);
	}
	else
	{
		printf ("No solution!\n");
	}

	solver_delete (s);

	if (st == l_True)
	{
		int pathway1[L][N], pathway2[L][N];
		int flex, fix1, fix2;

		printf ("Checking the frozen positions:\n");

		memcpy (pathway1, pathway, sizeof (int) * L * N);

		for (i = 1; i < L - 2; i++)
		{
			for (j = 0; j < N; j++)
			{
				pathway1[i][j] = -1;
			}
		}

		for (i = 0; i < N; i++)
		{
			printf ("%d ", pathway[0][i]);
		}
		printf ("\n");

		flex = 0;
		for (i = 1; i < L - 2; i++)
		{
			for (j = 0; j < N; j++)
			{
				memcpy (pathway2, pathway1, sizeof (int) * L * N);
				pathway2[i][j] = 1 - pathway[i][j];

				s = solver_new ();

				m = network_inferred1 (N, L, pathway2, S, g, r, s);
				m = network_inferred2 (N, L, R1, C, S, s, m);
				m = network_inferred3 (N, L, R1[N], C, M, R2, s, m);

				st = solver_simplify (s);
				if (st == true)
					st = solver_solve (s,0,0);
				if (st == true)
				{
					printf ("* ");
					flex ++;
				}
				else
				{
					printf ("%d ", pathway[i][j]);
				}

				solver_delete (s);
			}
			printf ("\n");
		}

		for (i = 0; i < N; i++)
		{
			printf ("%d ", pathway[L - 2][i]);
		}
		printf ("\n");

		for (i = 0; i < N; i++)
		{
			printf ("%d ", pathway[L - 1][i]);
		}
		printf ("\n");

		printf ("Checking the links:\n");

		veci lits;
		veci_new (&lits);
		lit *begin;
		int var, flag1, flag2, flag3;

		fix1 = fix2 = 0;
		for (i = 0; i < N; i++)
		{
			for (j = 0; j < N; j++)
			{
				/* check inhibition */
				s = solver_new ();

				m = network_inferred1 (N, L, pathway1, S, g, r, s);
				m = network_inferred2 (N, L, R1, C, S, s, m);
				m = network_inferred3 (N, L, R1[N], C, M, R2, s, m);

				veci_resize (&lits,0);
				var = r[i][j];
				veci_push (&lits, toLit (var));
				begin = veci_begin (&lits);
				solver_addclause (s, begin, begin + veci_size (&lits));

				/*
				veci_resize (&lits,0);
				var = g[i][j];
				veci_push (&lits, lit_neg (toLit (var)));
				begin = veci_begin (&lits);
				solver_addclause (s, begin, begin + veci_size (&lits));
				*/

				st = solver_simplify (s);
				if (st == true)
					st = solver_solve (s,0,0);
				solver_delete (s);

				if (st == true)
					flag1 = 1;
				else
					flag1 = 0;

				/* check none */
				s = solver_new ();

				m = network_inferred1 (N, L, pathway1, S, g, r, s);
				m = network_inferred2 (N, L, R1, C, S, s, m);
				m = network_inferred3 (N, L, R1[N], C, M, R2, s, m);

				if (i != j)
				{
					veci_resize (&lits,0);
					var = g[i][j];
					veci_push (&lits, lit_neg (toLit (var)));
					begin = veci_begin (&lits);
					solver_addclause (s, begin, begin + veci_size (&lits));
				}

				veci_resize (&lits,0);
				var = r[i][j];
				veci_push (&lits, lit_neg (toLit (var)));
				begin = veci_begin (&lits);
				solver_addclause (s, begin, begin + veci_size (&lits));

				st = solver_simplify (s);
				if (st == true)
					st = solver_solve (s,0,0);
				solver_delete (s);

				if (st == true)
					flag2 = 1;
				else
					flag2 = 0;

				/* check activation */
				if (i == j)
					flag3 = 0;
				else
				{
					s = solver_new ();
                                        
					m = network_inferred1 (N, L, pathway1, S, g, r, s);
					m = network_inferred2 (N, L, R1, C, S, s, m);
					m = network_inferred3 (N, L, R1[N], C, M, R2, s, m);
                                        
					veci_resize (&lits,0);
					var = g[i][j];
					veci_push (&lits, toLit (var));
					begin = veci_begin (&lits);
					solver_addclause (s, begin, begin + veci_size (&lits));
                                        
					veci_resize (&lits,0);
					var = r[i][j];
					veci_push (&lits, lit_neg (toLit (var)));
					begin = veci_begin (&lits);
					solver_addclause (s, begin, begin + veci_size (&lits));

					st = solver_simplify (s);
					if (st == true)
						st = solver_solve (s,0,0);
					solver_delete (s);
                                        
					if (st == true)
						flag3 = 1;
					else
						flag3 = 0;
				}

				if (flag1 == 1 && flag2 == 1 && flag3 == 1)
					printf ("* ");
				else if (flag1 == 1 && flag2 == 0 && flag3 == 0)
				{
					printf ("r ");
					fix1 ++;
				}
				else if (flag1 == 0 && flag2 == 1 && flag3 == 0)
				{
					printf ("n ");
					fix1 ++;
				}
				else if (flag1 == 0 && flag2 == 0 && flag3 == 1)
				{
					printf ("g ");
					fix1 ++;
				}
				else if (flag1 == 1 && flag2 == 1 && flag3 == 0)
				{
					printf ("R ");
					fix2 ++;
				}
				else if (flag1 == 0 && flag2 == 1 && flag3 == 1)
				{
					printf ("G ");
					fix2 ++;
				}
				else
					printf ("x ");
			}
			printf ("\n");
		}

		veci_delete (&lits);

		printf ("Flexible Pathway Positions: %3d , Fixed Network Links: %3d , Half-fixed Network Links: %3d\n", flex, fix1, fix2);
	}
}
